{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Contact of an elastic membrane with a rigid indenter: an augmented-Lagrangian approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual: 7.55796079282e-06\n",
      "Residual: 6.04663326501e-06\n",
      "Residual: 4.88782492198e-06\n",
      "Residual: 4.05065021525e-06\n",
      "Residual: 3.32603080899e-06\n",
      "Residual: 2.74757531832e-06\n",
      "Residual: 2.28568557151e-06\n",
      "Residual: 1.91675781292e-06\n"
     ]
    }
   ],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "parameters[\"form_compiler\"][\"representation\"] = 'quadrature'\n",
    "import warnings\n",
    "from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning\n",
    "warnings.simplefilter(\"once\", QuadratureRepresentationDeprecationWarning)\n",
    "\n",
    "ppos = lambda x: (x-abs(x))/2.\n",
    "Heav = lambda x: (1-x/abs(x))/2.\n",
    "\n",
    "N = 40\n",
    "mesh = UnitSquareMesh(N, N, \"crossed\")\n",
    "\n",
    "obstacle = Expression(\"y+pow(x[0]-0.5,2)+pow(x[1]-0.5, 2)\", y=-0.001, degree=2)\n",
    "\n",
    "alpha = 100.\n",
    "Nitermax = 10\n",
    "\n",
    "V = FunctionSpace(mesh, \"CG\", 2)\n",
    "#P = FunctionSpace(mesh, \"Quadrature\", quadrature_scheme=\"default\")\n",
    "Pe = FiniteElement(\"Quadrature\", mesh.ufl_cell(), degree=2, quad_scheme='default')\n",
    "P = FunctionSpace(mesh, Pe)\n",
    "\n",
    "u = Function(V)\n",
    "u_old = Function(V)\n",
    "du = TrialFunction(V)\n",
    "u_ = TestFunction(V)\n",
    "h = Function(V)\n",
    "\n",
    "p = Function(P)\n",
    "p_old = Function(P)\n",
    "gap = Function(P)\n",
    "\n",
    "corr = Constant(0.)\n",
    "rho = Function(P)\n",
    "dxm = dx(metadata={\"quadrature_degree\": 2, \"quadrature_scheme\": \"default\"})\n",
    "lagrangian = 0.5*inner(grad(u), grad(u))*dx + p*ppos(h-u)*dx + alpha/2.*ppos(h-u)**2*dx\n",
    "form = inner(grad(u), grad(u_))*dx - p*Heav(h-u)*u_*dxm\n",
    "J = derivative(form, u, du)\n",
    "\n",
    "def local_project(v, V, u=None):\n",
    "    dv = TrialFunction(V)\n",
    "    v_ = TestFunction(V)\n",
    "    a_proj = inner(dv, v_)*dxm\n",
    "    b_proj = inner(v, v_)*dxm\n",
    "    solver = LocalSolver(a_proj, b_proj)\n",
    "    solver.factorize()\n",
    "    if u is None:\n",
    "        u = Function(V)\n",
    "        solver.solve_local_rhs(u)\n",
    "        return u\n",
    "    else:\n",
    "        solver.solve_local_rhs(u)\n",
    "        return\n",
    "\n",
    "def border(x, on_boundary):\n",
    "    return on_boundary\n",
    "bc = DirichletBC(V, Constant(0), border)\n",
    "\n",
    "def solve_Uzawa(tol=2e-6):\n",
    "    t = 1.\n",
    "    for i in range(Nitermax):\n",
    "        u_old.assign(u)\n",
    "        p_old.assign(p)\n",
    "        told = t\n",
    "\n",
    "        local_project(p + alpha*ppos(h-u), P, p)\n",
    "        solve(form == 0, u, bc, J=J)\n",
    "        #gap.assign(project(h-u, P))\n",
    "        #p.vector()[:] = np.minimum(0*p_old.vector().get_local(), rho.vector().get_local() + \n",
    "        #                           alpha*gap.vector().get_local())\n",
    "        #direction = assemble((h-u)*(p-p_old)*dx)\n",
    "        #t = (1+sqrt(1+4*told**2))/2.\n",
    "        #if direction >= 0:\n",
    "        #    corr=(t-1.)/told\n",
    "        #else:\n",
    "        #    print \"    Restart...\"\n",
    "        #    corr = 0\n",
    "        #    t = 1.\n",
    "        #rho.vector()[:] = (1+corr)*p.vector().get_local() - corr*p_old.vector().get_local()\n",
    "        res = max(0*errornorm(u, u_old), assemble((p- p_old)**2*dxm))\n",
    "        print \"Residual:\", res\n",
    "        if res < tol:\n",
    "            break\n",
    "\n",
    "h.interpolate(obstacle)     \n",
    "solve_Uzawa()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAEACAYAAADyRL7nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXmQHMd95/vN6nvOBgacwUFcAxI8RclDciVZ1q6eBFpe\n0WtF2IC0lFaUQxKBcMjxQtSLAOXYPzZi40X4gS/W1IuQvYGhJVvSo/AoUM9BhSmLAihLlrQLWgR0\nAQRJEY0BQADEAHNhjr4r94+qmqmuriOzKuuayU9Ex0x3ZmX9uiv727/65S8zCaUUEolEIgkPJW4D\nJBKJZLUjhVYikUhCRgqtRCKRhIwUWolEIgkZKbQSiUQSMlJoJRKJJGSYhJYQcoyhzigh5CAhZI/+\ntxzcPIlEIkk/xC2PlhCyB8AogMOUUuLaECEnKaX36/+XATxNKd0n0liJRJJ+CCGjAPYCOAVgDMA4\npXSWt27UZYGglHo+tGqu5WMAjllem2FpWz7kQz7W1gPASdP/ZQBH/dSNuizIQ1SMdhSAVfWnCSFj\ngtqXSCSrAF0Tpo3nVPMW9/DWjbosKKKEdr2gdiQSyeqGxylzqxt1WSCyQRvQmYbmZptxFF9CyH4A\n+wGgVMzdv3XrRgBAPrcARWmiVl8HAMhk6ijmZ7FYHdGPo+gpTqJaXw9VzQEASoUptNolNFs9ehvz\nUEgbtYZmTjZTQyF3E4u1Yb0NFT3F66jWh6Cq2tvvKd5As9Wz3EYhfxMARb0xqLWRrSKfXcBidQSE\nqCCkjZ7iDSzVNoDSjN7GdTRafWi1SnobcwAI6o0BAEAuu4RcdglLtQ0AAEVpoVSYwlLtFlCq/d71\nFidRbw6g1S4CAIr5Wag0g0azf7mNbKaKan1Ib6OJUmEaS7VhGCH03tI11BpltNsFUKqgVJyCqubQ\naPbpn88iMkod1fp6/TNuoJifWf6MV9pYh3Y7r3/G02irBTSavaFdJ0JU9BRuCLlOS7Vb9DbCu05L\n9Q3Lz4NeJwAoFmZCu06/Pr1wg1J6CwLw0P9WolPTKlPdX/y6cQZAzfTSOKV0XP+fxylzqxt1WSBE\nCW0FNkZSSk/ZVdY/9HEAeOc7+umJH2T0kkFLzR79YWaz5fmw5blV73v1h0HGpo0Ry/N1lud9APow\nt7ANg30XHdrYaHlu/Tj69YebHZssz4d8tLHyXLPXOoY5YHle0h/2bWhYv6Pir9Pcwk4M9uVMr/i/\nTu5tiLlOK30BCHqdNDZYnou7Tj2bFy4gIFPTKn72fev57OnZPFGjlD7gUMzjlLnVjbosEL5DB3o6\nVxnoFlR95O44SzuGx5MGFqrWL3qySZO9abIVSJ+9CYLHKXOrG3VZIFyFlhAyRgg5qP9/SE/3MjgE\n4GOm548ZebTQ0iMeYzHAuFVKA5euvT9uE7hIk71pshVIn71JwcspY3Xgoi4LimvoQD/xKQBP2pTt\nc6gLUcZJJJJVyWO6A2fkqpqdskMAjkEPLXrUjbrMN6JitL7J5xbRHY9KJrcO/8+4TeAiTfamyVYg\nffYmCTenjMeBi7osCLELbUapx20CM4N9E3GbwEWa7E2TrUD67A1KGxRzanq+q0kj9kVljPSVNHCm\n8kjcJnCRJnvTZCuQPnsl8RK7RyuJlyi8lEGlEPo5JJIkE7vQZjINdOcHJpOBvktxm7AMi0AWeycS\ncbsnytYkCXaS+oIk+cQutMX8DNIitHfvOBr6OUQK487tR4S1FTYstrJ+NlEIchR9QbJ6iD1Gu5ii\nxO8TZx4X0s6cWnd8iOT02YNC2wsTkbZG8fmK6gtpoU2BOZXtIekmdo82VbgvydtF7LftNPbfUXYi\nstXpmnB7wZx9QbK2kULLA3FeJB1IgLBaISlyL2K21XrtPIXXoy9IJGZcd1iIgrF3FijrYhVJI3HC\nKgmNJA3E8dKzeeKkyyIvTNx7X55+5wXrwjf23LntauDzrTZiv7esNawrMCWXVyf2hRZPDYPzF9KT\n65l0W63X/dUJuUuThJ3YQwfGWppJxvhyzSzcii0x28LD4uJ21/IoBy4GPX7SvWxNEnNqHTMLt2JO\nrafa05VER+xCm1TS4LHaYRbPNpIzCuxlh9lWL1FOEkY/We2C2wLBlFqM24zUErvQlgrT6F6wOD7c\nBHZ0xzcitMQZVvHcuD0Z9rJgttXt/SVFhK19Ya0IrsQfsXfbtpqMjskSd11YHI3ImhWC5CnWYrDX\nL6y2JiVv06kvpCV+L4mW2IXW2N8oTli/GJPXfy9kO8SKyOyNcO0VSRBb4xBer74gBVdiJvbQQZwk\n4YuQlBjqasL6mcYZbpAhBQmQAKHN5xbQvdlf+PgR2ZHhHwk4b+AmmFl3y4+iO1lAwrTV/JmLEl3e\nvpD2DIU2VTCrWjfglLASu9AqSjPS8wXxYkvFqz7P6fuUgcj7tDcOorJVlOj66QtpF1uJf2IX2lo9\nugkLQUMFExcfwb13/18c5wt0usBcu/QIdtzZaW8SUnSGlFrXa3a2hk0Q0eXtCyvnlKEEVvTNEfdi\nZf+ucUrpLG/dkMoOATgIYBbAKwAOUEorTu8ldqGNiqjisXGLK7Aipi2qJEJYrdjZZNhqJ8JREHUO\nr/RumThKKb0fAAghrwB4GoDTlDy3umGUnaOUfWWh2IU2k6kDCDf2I0pk+/vedDmHkFNw4yakuT7H\nH9jEYdhq936iFF9WwXXrC+znkmLrBCFkDMC08ZxSOksI2cNbN4wyP8QutMX8LMIUWpGe7Lat/79N\n+8KaZ4LHQ+3d/N0QLRGLm63W9xyF8HoJrl1f8HeeVSm2G3QP0GCcUjruWNueUWi35WamCSFj+k61\nTHXDKNPPXyaE7NXrPATgL53CGkAChDbMhb9FhwvOnD3YEZeLQmSD3PrPvvEFrLvjrwRaEx48tkYp\nvE6Ca+0Lq50WMphq9bFWvyFg9S6eXVvd6oZRBnTGa6cBvATgfqfKsQttGglbYJMYV00y5s8rLNEN\nM4a7Sr1aWwgh+wHscqlyjFJ6HNpte9lS5iR+bnXDKIPZe6WUniKEjBFCyk5ebexCS0JaQDmMwS9V\nqYcmsmGIK1Hin5DBiihbjc8xTMEdVABF8Ge7VsSWI4RQgY2w2oQNXOsSQhBC2RiAp42BMlNZckMH\nPcVJAGIX/g5DZOdUYPvup4S3G6b3Wr79r13Lo0xALytL7uUetvISppc7pwJbQugLa0VsWTCJHYDl\nVKvjlufTlNJZt7phlEET9sOmsj0AnnN7P7ELbbXOE4rxJiyRBYArE49is6AVvKIID9y88AjUrc+H\nfh4WnETdEOCbFx7BQEi79obh5V6ZeBTY8Y3ErCa2SnmMEHIQK3msj5nKDgE4BmCcoa7QMj0DoaKH\nQQAtFGI+rovYhVZVc3Gb4Ig1TNCoBfe8wxZYs6C1alviv8AeGPa2aluW//fyfv0iUnCNvmCEEkSR\nVK+2TRXMtaOdgquHCYxQwXFL2T6OumGUdTz3IunfQy5EerOiY7FhCexqnH9ufU+ihVe0hytabCWr\nj9iFtlSYAjActxkdOIns5h1f5W4rDIFlFdfMtmeFnzss3Gw1v1+RohtEcK19QWRWQlK9Wol/Yv8d\nbrVLQtoR5c26ebILN+9lbmdKLQoV2Vm1Z/nBijq/W9j5w4bVVj+fgxd+rpVTX0jCFGxJ8ohdaJut\n5Nz6en1Jbk6/m6kdUQIbVFTozO8EOv9Uq4/5ERQ/toYhuKy49QURYpuEtZIl4og9dCACEZ1SxJdD\npMCGjQhx5GlvKLsg9HxmRIYWws7DlaxNPIWWc6myMQDG1LsygOfclg4DgHxuHt0TMKKFVWTXDx9z\nLBMhsqIFVrnlJwDEi6ofnGwwBNiwNSiiMhe8VhJz6wuAmAGyJMVq48g6WE2weLQ8S5XtoZQ+aTwh\nhBwGcMCtcYW0GU21J8pbrGxuzvb1oCIrWmANUcuQFtoJEFk3rLaK8nxFCK6b2Dr1BTMyG0Fi4NoN\n7JYKA+C2VNgBQgiXe1prpMObBYDJy3s7ngcd8BIZY7SLlxavfUhI21Fg2Coy7gsE/xFzur7WvuBE\n0JCUjNWuDrw8Wp6lygBttsZ5QsgT+vMnbOokhiBfgiR4sUkICYSJ+f0F8XSDerdB47bSs5V4CS3X\n/FhK6bju0RrhguPoFmpjBZ/9ADAyUsaJ0/8VALB15CfoK13D2QnNWyj3V7B723fxr2e+AADIZBp4\n8K6v4HTlE1hY2og2KG4b/Spm5+7FjSltFHjTxmPI5eZw8ZLWxsDAWWzZ9H2cff1x7Q1nF3Dn7q/g\n1XOfRbN+CwBgy+hh3Jx+APOz2hoRQxu/B0Vp4vqVjwIA+gZ/g/KGH6PZKGPitS9Bzd7E4OjXMFf5\nDNSm5pEP7jqM6o3fQ2PuHgBA76YXQGkWS29/GABQKP8ShXWnMHP+c9qHkJ9CdvuzaFU+DbS1Ldcz\no38LdfIDoAu3AQCUTd8DWv1Qr79f+9zW/QI3S2+hdFmzq1i4jtqWF9Az8UkQfYbd4s6vo3Dtg8gu\nbQVp9kOpbkKmvgH5ae29Ndb9Eu3SZZSuPAwAaJeuorbpRfRWPg2AAKBYHP06ilc/jEx1EwCguvkF\nZKpbkJ95l9bG+pNoF26gdFV7b62eS6iP/BC95z+t9QOliaUdz6B4+WFk9M+4uuV5ZBd2ITenpUXV\nh14GzS4se7IUbUDNoXfik9rzTBVL259F6a2PQmmsQw3A0tbvYHD+NtC5d2ifz/APAaUF9e3f1z6f\ngdegDJ1A+/yfap9xbg7ZHc+gNfFJoDmIGwAyO/8e+en7Pa/TzfOf0a5J4ToGdnwTs+f2Y6bVhyxR\nse32pwC0MfHalwAAw1ueQ6s5iOnJh7Q+t/5l9A2cxpWJz2qXungFg6PfwKuvPQ5V1WKu99z1JC5e\n+mPM69d6x7YjqNY24drkB7Q2b/kp+norqEw8igwIBvou4e4dR3HizOMAJQCheM89T+HViX24ubBV\na3P0COYWduCtyfcCWPk+AROQxAuh1Hn1LH1h2wOU0odMr80A+JCdR0sIOWjEaHUxfYJS6rYkGsbu\nK9KfvbjJl/FBbqv8eLNqu4AZwrx7RRdBvFhf3quaAxg2vwxzkGMww+hFMtpqJoiXGyR2O6TUoLYL\nUDJ8/S+IVxtkUKxn88TJoOvDbrxnPf1P33rIuyKA//aubwc+32rD69IzL1Wmr2BzylRnHMBzepzX\nkcVa9LPC/IYMKm/8H76OCxKLDRKvNLxDM3Ptnq5HmLCez85WL4J8NoGuiVrExd8+zn2cnMywdnEN\nHfAsVQZt0GwP7BdmEI5fb9ZvZ/cbkw0isCJIYkqO1SZmr9cB47Py4+HOqj2+vNsWlUFXCTss6V1M\nS5XpojxqWjqsDMBzsj0hKoAMn9URY4gs4fwi+xFZEQJrCFkuJQt/z7V7kFPqmGv3BBLdKZ/pYX7E\nlrcvGPgdGEtSTq2EH0+h5VyqzHXxWzt6itcheuFvJ/x4s2ZPtryLbXH4OLxYO8/1RkLWomXBsNX8\nPvyIrl/vlldsy7vGY90eXZIuYr//qdaH4jaBmZsTn/Ks49eL9SuybnHPoSt/4KtNu/bdHiKwszVI\n+34+T55rZ/QFPyGlNMZq21TBzWaR6SHpJva1DlSV3wQ/8dmg3iwAtPVUJSeiChWwik+WYTKICKF0\na4PVK3Wz1Wif18P1492yerbmviA923DgnP7vWNdvmantY+bMK17bgAQIbVLh9VR4RTZMgQ27jaDn\n8xuHDSK4YYhtxzk4xdZPrHYNxml5pv+71fVVpmdSjcJ+NiyPbfELbU/xBoCRUM/B6806iezAzq/Z\nvp5Ukb2x5R8Tn3VgiOaNLf/IfTyP4PJ6t15i69QXJGKwm/6vCx9XXb9l+vPjer3ljRh5bTOIPUab\npPVovajPdKcEhy2yfuOUc+0eKHPsC5XHhfH+/Njq57Ph+fzdrq1dX+C9C0pjrDZCHKf/c9b1WybK\nNgApFFre+KwobxYA6rPv4mvM2jbHl9yPiFgHqAYXRrmOj5PBhdFAPyo8iBBbp74Qxe7GKWADIeQV\n02O/9yFd8Ez/d6vrt8zv+WyJPXSQZni8WV6R5UFUeGCm1SukHQBYl130fayfFC/ecAJP3NbvpAYW\neGO1ccVp25Rgrsm87dQNpym4uui6Tcs/pt+yT6N7oWongXOr67fMDe7jYhfaQv4mgHVxmwHA2xvp\n2fji8v9hiKxogZ1c7zwpT6SosrbvJr5OtvIKKM+khyBia+4LXe3KLARH9Kn5LDBP/3erq89s5S4T\naBuABAgt4LyoTVB4wgYst3yEtADEL7KsdallUfWwxdUL8/mtomu11QqP4PLU9Su2Rl+QhAPP9H+3\nun7LgthmR+wx2npjMG4TmFm8+nBqRBYARqYexEyrd/mRJKx2jUw9yHRcGJ+VnyyQxasPu7fJEauV\ng2KOPEYIOaiP6O9F9/T/jzHW9VVGCBnTlx8AIeSQJbPArc0uEuDRshPWavNhDGCIFlkegTHEKy0L\nn8y0erls5fVuRXq2YcZrJZ1wTv93qxu07ElYcDvOjtiFNputAhC/U0AYXoI68AbT8jdxiazVa53u\nucx0nIHI6ZMDOb4Y5XTP5WX7WQfSWEU0DLHND57xbo8jVpv0XRjaVMFCMx+3GakldqHNZxcQhtCy\nwurNzqo9UIZOeLcnUGT9CqzB1f7fOh4T9px0a/tewmu2lUdwWb1b0WJb2vBTzzphsQZniKWe2H9D\nl2ru6wckieUtUhxIksgCwD3XPtDxPM6FP7zObbUVAFdsWeRnysL0uT9jqifzaiVAAoQ2DFjDBjze\nrChECQKPCCVtVSVewWd9r6I+W9YfTNHbxEtWL7GHDghpI+kLfy+Tm3MsYvlyihRZL242i1hU6r4F\nliM5fZnBXJX7GMO+OsNt/Uyr1zOcwBIiYKnjGUJw6QtdbTHGapMep5X4J3ah1RaViWbhbz+YvZbs\njmds60QlsjweLACcZIgpA/5ElbUdVvF9ef3PgWbRM5Ybpdi6YfQFmYUgYSH238+l2gameqypXaLD\nBmZaPjYQBKITWevt+P1T73E+X7O0/AgT1vMYtrKEFFhCCSI+c7cfUN6+kPZYrUoJFpt5poekm9iF\nltLkhg26YnDN7skVXt5sFCLrJE7Fdqe4RSWuTrid32orq+C6nk/A4Jfj9TX1BaExfDl5YVUSe+gg\nzYjcSNEJFpH1PIcPYeXNmezLNbjqGzZ5hRZueoQTvEIJXiGCoCGEuJApXukido9W25wxWlhu4+y8\nlMzOv+c6j5eIBhFZFo/vpcFTzCK70Mx3PHixHs/ahuHh/nzoZ451wvZs/YQQrH2BxatNe/hA4p/Y\nhbYhwCs0CPu2SzXFPIOGDIKKrGvbunjdXt3mWi+IsLLA0/7I/B2uPwpePyxhhxGs11t1iX9LJFZi\nF9pWK9p4oV9vFgDozTuF2BCWyFrjn1vqw111whZXJ7zOa9jq5YGHJba8QiyqL6QFVSVYbOSZHpJu\nZIzWB0G82TBF1g0eYfU7ctzLGKc1bHGK63rFb93itizpX054xWu9cmtZUr1YcmplPu3qI3ahLeTn\n4H9HCfG4xdqUjT/wPD7ILapokf1V7xtMAisqJceuHTfxNQvur3rf6Cqfa5aEi62owS+WviCRGCTg\nd5N41ghreURu1GygTAM3EfYjsm6pWgvNPJotdwGNIu+R5RwLzTwyDl3RK27rhNvnGeSOY/n6q/Y+\nipyWK7EjdqGtNwYiO1fQUV918oOu5X6/wH5F1gnDUxyr7ewq400sFxWT8zrv3Qu3O3rffsU2TLz6\ngkRiJnahFYWIjAMvbySMhbTDElkrLOLqW0Q5j3Gzw01snd630+cUllcbFJnmtfaIXWhz2fQkizcH\nzjqWif7i8oqs3Yh+JX8NgLuwhTFazNKmVfgNWwH3gbskiO1Uqw9k8DeO5SLCByyOQ5QhNZUS1Jo5\npoekGym0HDQHXxXanpMI+BFZO35Nph1FllVcg36pvETXENw3TUIL2P9w+CWM/dKU8q+FtynphBAy\nauzLpf+1bvHNVNdvmanOMZvXDhFCKCFkhhByTN+g0ZHYsw5YF5UJitftmpcXMtXqQ++lP8Hi6N93\nlUVxG8orsovNPD5Ruxvf6vll5+se4srrkdjVL+aa9jbp5+7Nd2ci7Jl/J54f/Neu1xea+a40MKds\nBK/punb4zUKoT3wKxd1f4T5OwsVRSun9AEAIeQXA0wD2+ajrq0zfeHEUgHlTRoNzlFLvkXyd2D3a\ntQqvN2uHm8h2vebiVYq+7fNqz827tcPufYoMITgR5EfS84dbxmldIYSMAZg2nlNKZ2EveK51/Zbp\nz49TSsdFvJ/YhVZRWnGb4ImR0qPmZ7rK/HizIkIGduJjjXvO6onxXgLrRr2R9Xy44XQOq/Av2xqi\n2DrhR1Dt+oJEKKMAZi2vTeviyFPXb5kXZULIXj3kcMgtrAEkIHRQKkwh6MLfUS0tV731+WhOZIJH\nZK0cVSqAjci6iauXcLIcU8h3/3ga57SGFRYbefTmG/he8bWV15p55llmbpMarPiZNeYUWqje+jyq\njBs5rgYoJTx9Y4N+G24w7sMz5JnF5FbXb5kX47oHDELINICXANzvVNnzk9ODvHuh7WE+Zj6BQ/29\n5ueU0ufc2k/C5oyso8Q9Fz6Ope3PMtUV4c0GEdnFRh6PNO/Akdzry685CawfcXXD3J5VdGvNnK3Y\n/qf27fiH0soW3nZiaxevdYI3Xssbq+XpC2uQG5TSB+wKCCH7AexyOfYYpfQ4tFt6q5foJIxudf2W\nuWLWQErpKULIGCGk7KSNLN8w5oA0IeQggAql9DndlX4JgKvQ0hByU60EiYeZZ4IR60Lagga6RCfd\nG7fkJdPltRNZN4FtNfgXZM/m247nMAuunXebV/PL3q0Bq9iG7dXaYe0LdgTd5mY1rnnA4dlWYCN6\nlNJTPHUJIfBT5maYHlp42tBF03GODqjrZeQMSJcB/IXhwVJKZ62GrBV4vVnbNnx6s3bxWKvIOsVW\nW43M8sMPbsfbnc8pdtvx3OY9s8ZrRf2AeeXU+kUOiDljFTv9zvq4+bkRF3Wr67fMgwqAw6bj9sDD\nofTyaB2DxTaq/wCAih46mIUWZniOUlpxO0FvcRLAJg8zksGiw+aMQRAZMrCK1DezZ21F1oqTsKpN\nNsFVct2erNGm2ct18m6LuSa+mXWeDMITs2XByavlCR+E0RckXTym3yUbYcvHTGWHABwDMM5Q11eZ\n7mgaGQqHoIc1KKWzhJCKHgYBtFCIuc0uvISWJ1g8qhtqGPIKgJOwicfoBu4HgJGRMk6c/iIAYOvI\nT9BXuoazE1qYt9xfwe5t38XpV78EAFCUOu6+8ymcO/8oqlVtAO220a9ievZe3Jx+t2bw8DFkc3OY\nvKy10dN/FsrwP2P2zT/Xzp1dQHnXOG5OfArt+i1oQ0Fm+zNQZ+8DnXuHdp7hHwJKC+rbv49eqqDZ\n/yaa606i99znQHPzUHPzeHvzP2Ho8h8i29Q8muu3Po++2ftQWtiJDVBwbejnIDSD4WltAPPtnrew\n2HsBd0++HwBQyy7gteGf4cEb70Ne1QTyxIZ/wW3zd2JdbSMA4FTfWZTUAu5aGoVKCd7IX8Wl3BQ+\ntHgvVEpwQ1nCD4pvYN/SfVD0vdf+LnsGe9rbsJX2o58W8P/iHIZRxHtwC1SV4GUyjQt0CR/L3ApQ\ngotqFc/hKr5Q3AkFBCoo/mr+IvaVRrCtoP0IfGvpbezIFvG7eS2c9S/1GVxTG9hXGgEAVFpVPF+b\nxON920EIUIeKv65N4JH8ZmxS9DbUC7hD6cf9ZB3QBn6mXMMCmvj32AI0AAoVz+Rew6dadwFNoKE0\n8A+lM/hI7U6Udc/veP+vcFtjBKMN7byv9v0Wbah45+JuAMDlwiSu9b+OB6fep33GmSpODp3Au6cf\nREH3Ss+M/Aib5m/HcFX7cbdep3r/61gaeB0bLv8hAKCVn8XU5u+j58LHl8MFizueQc/EfwSB9pmr\nm/8RaPVDva5dW7LuF1D630D74scxAxWZ4lUMbD+C2d9+HlTVtp8p7/4yFq/8EeYXdwAARrYeQaO2\nCTPXP6CVb/gpir0VXLrwKACgt/cCdm4/gtNnDwJUAYiKe+96Eq9O7MPNha0AgHtGj2BuYQfemnxv\nx/cJmLB+Bbmh1F84Kdg56SloAghYPE1K6T6OukHLnrSxjcXzXYZQSp0LNe/0AKX0IdNrMwA+ZON2\n7wFwmFK6y/QaBbDLzau9Y/c2+ssfO19AlmmGXlkHbrdobgNh1tvC3sqfLk9YcLqd5Akb+PVmnQa/\nzNSaOfwZvQP/nWiDYWZP1u4L4+i91hmChAX7C2D1dK0xXLNn+2f0Dvxd/nRHuXVig51Xa43X2sVq\nnQbF7LxaN4/WXGbuC26ZB14xWq+1ab1itCz7hvVsnjjpNDjFSmF0C93yf36eqe75T/7nwOdbbXh9\ni3gD0lYcg8MScdiJrBk3kVWbmW6RrSsrDxYc6lvbtZ7bGsaw2s0Sr7USNFbrZ4DTLU4rl02UAB5C\nyxmQrgCYNZ7rfyteMdpiPlwtFjngUBt5yddxQebae3mzbiL7T7jsKbId2IglaSieD7c2rEJuHSwz\n7PsnXO6yn4Wot+Qx8NsXJGsTlvQunoD0PgB/QQg5By026zQveRmVBov7RDVZAQCI7rnwhg3sYA0b\n8GAVqWKrAEC7NXUVWRtx5cFcn+bVzjb1sILazHSEElqNzHIood7Ioi/nPF3XK+XLCmu6V5BULyJw\nU1HJ6sdTaDkD0hUAT/AY0Gj281SPlcLUu9EadB4dFw2vN2um3sji32Vuwa/acx0ia+vF6lgFVmmw\nrZmh5lfi/KShrIit0b6D2Jp5nzqC0xltWqt1UoOX2LJMZOCZwOCUfWB+XVRfYNlDTJJ+Yp+CGyc8\nA2F+YQ0biPRmvQa+tErOXqyTwJpfN4ur9XWjrQ7v1kZszV6t3ftxWgVMEgOUMKf7SbqJfd5Jutaj\nPe1diQGWwZkg3qzBSdq58MnyF8XixRrCqDRIh2gaz62vu5WZn3d4yJa4rYHxQ3CSzrjOVONdmJz1\nh8tv/NyjnPqvAAAgAElEQVTaF0T9MFuJMjQmCY/YhTabYZs6mQRafeeExGdF4+TNnqmv/IjZeSNO\nXqxZLJWG98N6jLm9jgEzF7F9XZ3vst9tYMz6w8MyKCZyqnOr7xxz3TAzDxKzcanEldiFtlofitsE\nZkqXPxpKu7xhA1bv7tHCrQDsB77sRNZOYA0yje7HyvH2gmsbgnBIGftEZhvTe4oKrx/OsPqCZHUS\nu9CuZvzelrJMUDDgis3aiJxVFA3RtBNVqxdrV+7Url0YoeMHwLRYvZNXy5tXGzTubdtmjHcukvQS\n+2CYojQBpCPI3i5cD9xG0NtXFm/WENmras0xLgvYiyxgL5wGmQbQzneWqflOD9d4TWkQqHm6/Lcr\nIwErg2NX1eAj7zzLKAZFRF9IFSrYJ7BIuoj9kysVpr0r+UT06ki1LS8IbS8oTt6swZHGla7XeEXW\n7LXmFuy9XauXaz7e0bO1fGmPNK5EMpfe7ofOz51H0vqCJNnELrRLteG4TejCaQS5Z+KTws/ldXvL\nMu3Uic8Xdmr/OHgidiJrDR0oHa9RZBoUuQVqeX2lrtGOuV0mW4s7Op6zhg/iIoy+IFm9xC60HBtJ\nxg5R7UfBRcbtgk4pNXuFeXR+tlZvVvtf+2u+9c/YiGt+niJTB/I39b/zdLnMKtTW9ry8WrWZQcFn\nVwzyQxQEp74gkdgRu9DGRRoX+xDpzZnF1uyBZiwCm6lDf1AUbraRravI1Kn+MI63F1ujTR6Chg+s\nP1RhDIg5EVYurST9xC60vaVrcZvAzJVtRwMdLzKP0ys+qzYzeGrhQlc6l503a/xvHdAyBDZTp8jW\nVSj6I3+ztSy4hodrFtvl4228Wif+av6Ca3kcuN2pLO78eoSWSNJO7EJba7ju0psoyvrCzmnho0X+\n+PeyN6uLbFYX1txcE5laG/m5BjK1dpfgmsXWzat1Ch/Y2ZrkOG3h2gfjNiFaKGFayY13QaK1Quzp\nXe2298LFSaGwxL4tepClEUUxmrW/bbYLG9hheLGZWhuZehukpi3Una21kKln0S7ot/mF4F+u0WwJ\nsF/2IJFkl7ZCzsmSsBK70ErE4Te+2RWb1b1ZAMsim5ntXJNChZb93C5mkK2r0G6OCNoWZ9PIu5U4\nI1fwskdf/3ovVpZoHXfaadatboCyMWh7IZYBPAjgCWN9bR7bgAQIbbEwA2BD3GYwMTPyz3GbwMXR\nqr/4tzkmC2DZk0WtDlRrQKkIBbrYWrzaTIPCEFyloU1cMDAmLjjamqK7zuqmF+M2YS1w1NhJW9+D\n8Gk4r3HtVpe7TN+44AFje3R9q65jWNkDkce2+Lu2mqI0mWyDZ6/K4ARNXRpROo93GggzkzHdDxve\nLIBlkaXVmia2y/W1ciNWu/w6Z7aB2dawJy6IGJTM1NPhHKQV3Ztcns2ke4t7eOv6LYO22ax5be1X\nAIwSQso8thnELrSNZnpSYvpn3hnp+YJusf1vC+u4jzGHzNtFk+AVtQJSKgKloql+xvT/ipDzhgvM\ntjqtUZsk8tP3x21CtNDupTGdHgA2EEJeMT32ezVvwyi69xyc1kWOp66vMn3Dg4dMrz8AYFYXVR7b\nACQgdCCJDmPdAe1/Z69WE0wFWWhCmjVCB+sGNc9WF11a1LqPWlDQMg2IGeEC1SK2TmGDqGHdaUHi\nmxsCdsHluX10q+u3DJb9Dg9gZRsv7lvb2IU2n1sEMBC3GUwslMUs/B0V/6PBtvGlIYh2wtsuZpCp\nZ7G8FExRq2yIbLuYQaugoF0gaBeAdt7fTD9WW5NCY90v4zYhleje7S6XKscopceh3Zpbcz+dBM6t\nrt8yq83PUkqf4znOTOxCm1HSkyTTKF1lrrsuuxh5ilc23+6Ib0607D03s2fb+ToAEABU/6sgX1fR\nGMwjU8isxGuhebrtYgaqR2oXawhholVLyyJuAIB26bKwttZSxoExuMRABTbiZd2Z26suIQR+yoz/\n9UGwii7+fmwDkIAYbbUe7QBTENZffci7UoL4RM9GzzpOQtguaCGExkAWakFBu5hBYzCPdkH72xzM\nLYcMzN6sNVzgib6XmJetSds/rHTl4bhNWNVYRUtPpzpufq5nBrjW9VumPx8DMG2ILCFkL8txdsTu\n0a4lBnI1odNwmSioQF1bB5Y0FMc4badAEmQaVB8Y0zxbFBRk6+qyB2sILNA5gAZo4m1ur/N/LU5r\nXZeW2EQcCvmW7Vsy74gbJnY74Uoi5TFCyEGs5Ko+Zio7BC3dapyhLneZLp4n9f+NuhUAzzG02UXs\nQpvJNABEt/CHQVlZ4l9Ypkfc7aIfevON5emnxVxzeVpqId/qWu9AybVxse18S2oW3LZp4e4VD9cq\ntoD5BsjwYrVjVoTbLLJtG4F14qKq7R0XNOMgsoW/bcJIQ9mFSM4dB4Sybz8vCt1zNLzH45ayfRx1\nucv0gTDHN+zWph2xhw6K+RnvSgmh5pCkLtLzCSoUZqFimbBgFkXr/+080R4FzWttDBA9nLASKrCK\nrEHbJvOgy5strHi1zzXY499mvFLgBnPhbP7p1BckEjtiF9rF6kjcJjDTW/m08Da9hCBILu0X+/UN\nD3VBMwTOELyVvyvHWMW2Q3DzBI1+0iWwZk/WLLgraV7u3qySa+MLxZ0dr5nDBkmLzwLh9AXJ6iV2\noQ0Tv6O5zreAyVqk3CxAdvFMxcZeJ7EFNIG0CmY7DzT7Vl4zP5p9nccYbWjtwvY8dt6sYStL2CCq\n+Kw3yeoLkmQTe4w2XcSfcG+O0zphpHmpZnv1QTEnzLm0RszWK4PAWu4lsnYoOU1c1Yg+W7vJCuuy\niz5air8vSNJD7B5tmhb+XhzlW+zZ7gvMMitJ1IDOl2vnl4XMjJ1Xaw0jtPPWwazuh4G5rlu4wMmb\nBYCvqG8u/88aNrCGVayfW1jxWYC/L6Qe2rl/nNtD0k3sQltr8M/Hj4vi1Q+H0m6QOK1T+CCbb2Nv\nfhOAFa/RLHBWsTX/r+a7BdftYT3GTsDtRNb8I/DHyhbXzwCINmzgNMBpvB5WX5CsTmIX2naKFivN\nVDd5fgGjwE1wzGK7LWOTs+sgtk7erZ0H61Rmbcd6HjNmkc3m27iVlLrs5xkEiyqtyyBT3RTp+STp\nJnahldhjFQ5Wr9aJjhCCjdgC3d6t9fbfTnSNenZCbW3fOK9dOMML64+Ln2wMcfHZbtxyaMuK/x/h\nQY9v6KCSnh1K1jKxC22pMO1dKSFUN7/AfQxrnJY3nsji1X67/VbHSL6b2Fq9W6twOj0M7AS2I1xg\nE5cFVvJ+v91+S2hKV5jxWcBfX5CsXTyFVp9TfJAQskf/y7SbIiHkMEu9thrfL7Kbp2HnoWSq3nHE\nMOH1arcTbeYbi9gC3bf4TqLqVtYhsDbnsIYMDHblnPuBlzcrOmzAEgaKuy9EDVG1TBSWh6QbFo/2\nKKX0SX1hhXFoWza4oq94w7TYb6MZbIUrr1srkeRn3uVaLjpO6yUgVgGyDoy9W1lZsMdVbG28W7u4\nqpPomo/rwDLw5SSyhXwLD5q2MwrDmxWxBq35+nr1BVbW0spdaxlXmfKzZYPu8U6jewXyVQOvoIoM\nH1i9OTexVZROQbSKrZvgAp2iaxZSp9ed2rLGZN0mJlhFNkxvVlR8ViLxwssf5N6yAcAet3UZreRz\n4S7EIdJjaKw/KawtVoIIyQlc75oxZhW5roEpQyQdYqp2nq7TcV1ibnN+w74TuM4tsnaE5c1aMfeF\n1byYjEQMXkLLtVisHjLwXMmmwwDF/TYxzlFV6xeoXbghtP2wvdpJaO3bia2rd2tgFk+vhwkngbWG\nC8x2zfkQqzC8Wda7Fda+ECTjQLJ68JqCy7xlg75+47Tb3uamuvuhx3CHh4dw4vQXAQBbR36CvtI1\nnJ3YCwAo91ewe9t3cfrVLwEAFKWOu+98CufOP4pqdTMA4LbRr2J69l7cnH63ZtzwMWRzc5i8rLXR\n038WyvA/Y/bNP9fOnV1Aedc4bk58Cu36LZpB249Anb0PdO4d2nmGfwgoLahv/z56qYJm/5torjuJ\n3olPQC1MQc3NA1u/g+zFjyGrby55/dbn0Td7H0YWdqINBdeGfg5CMxie1pz/ub4KLpSu4u7J9wMA\natkFvDb8Mzx4433I6/lSJzb8C26bvxPratoi2Kf6zqKkFnDX0ihUSvBG/iou5abwocV7AQCTpIof\nFN/AvqX7kIOCNiX4u+wZ7Glvw1b0YwNK+Cr9LYZRxHtyt0BVCV5Wp3GBLuFjmVuBEnCxXcNzjav4\nQnEnlCKBCoov187jT7Kbl/Nwv7X0NnZki/jdvNYV/qU+g2tqA/tK2oJAlVYVz9cm8cX+7QCAOlT8\ndW0Cj+Q3Y5NSBAjFt9oXcYfSj/vJOigKxU/oJBbQxL/HFhBC0d/K429yv8KnWncBABpKA/+AM/hI\n7U6U1SIUQvGDzK9xW2MEow3tvK/2/RZtqHjn4m4AwFTpMi4oFTw49T7tM85UcXLoBO669n4U2trA\n4JmRH2HT/O0Y1vNgrdep3v86lgZex4bLfwgAaOVnMbX5+9j01n8AaWu5vos7nkHPxb2guXkAgLr5\nH4FWP9Tr2rUl634Bpf8NtC9+HDNQkSlexcD2I5j97edB9cHf8u4vY/HKH2F+cQcAYGTrETRqmzBz\n/QNa+YafothbwaULjwIAensvYOf2Izh99iBAFYCouPeuJ/HqxD7cXNgKALhn9AjmFnbgrcn3dnyf\ngAlI4oVQ6jxnWw8RPG3sX66/NkMp7ZrOpa8+bhbhw9A2NDtu2eSsgzt2b6O//LH7HiZzqvt2N3MO\nd7MGU6r7Yttu69JOtVZ26e2t/CkWR/9+5bxt++OcXnfa2sZuMfC5ZvcavQuW7cfttiM3r4PwmeY9\n+Bu80VXHunYtEM4W33axWKt3bXjfn2neg6/lzgBgCxnYebM8YQNej9b6urkvBMmh9Qpticij7dk8\ncTLoZomljVvprke/yFT3zP/9xcDnW224erSmfXUA2G8nAd2LNW1cZpQdZtkfKJOpA+BcgDsmWj2X\nmOoNZpZsxZZnH7HBXNVWbM305hpdYmtedOYSmUcxqwmZsUg4sCJ2ZsE1i2IQ0XUa6LJbXcwck71E\nNO+QZZpt1CJrh9EXgsRnZcbB2oElOeoxI48WwF50byfxMXNlQkhZ3+IB+nGjbo0X88lOTjB/keoj\nP+woE5XOxTpYYycwdt6eIVbHMxeXX7NLmXLaKsaIp1rjqn7qWmOxTvYcz1y0Fdkg6/GKwu46W/uC\nHTI+GwyeHH63ugHKxggh+/XXj5q1jBByiBBCCSEzhJBjXjrnuUwiz3YS+muzAJ7UH55EsfD3kFLz\nDB+w0Hv+0x2hAz8E9Wr7co2uEIKTZ/vxxd9Zvh0HOre/MTCLoF1IAfC3vYyTiDvlyH6ufTe+hc4t\nvJMQMnBCRF/wIsoc8YRy1AhbEkJegZbD36U5DHW5y3TBfcC4K9cdzWNY2Sr9HKWUeVFieSkRjufB\n+8XlSUFi9WwzpDv+Xsw1HcXO8D6dRNILr+PtztubbzB7sqwiKxK5QWM88OTwu9X1WwYttfUJ02le\nATBKGGfGWol94W9iIwa8DCreA2JBGMouYKrVB2qTiuYUj3XDyau12yXXKVbL4tk2oS6LmHWxcEP0\nrB6ugV+xtcNJ2M0C28TKBQwqsmF7swBAleaayp8llGt67QbdOzQYZxmvseCYw2+Tp++W7++rTB+f\nesj0+gMAZk1ZVWU9AWAWwEMA/tIt4yp2oe0pTgLYHLcZTCzteIarvh8Rtm3Hp9ge7fn1yusOOzOY\nRdBJdP3gumC3jQd7tOfXjvFYnnxZXpH1y9KOZ1z3bg6abZBybgjIOuDJ4Xer67cMlmypA+gcnxo3\nhJUQMg3gJQD3w4HYhbZa55oTERos248XLz+M2pbuVZvC9moBfrEFgPfN34sfFFfSu5y8W4Mgosu6\nPoFTVsFHGrvw49zZrtedRFZUyMDNm3Utu/phYOt3hNiwltBz6He5VDmmr6vCnMPvUddvmdXmZ82Z\nVWbvVfd+xwghZSevNnahVVVxXpQbQQfEhrILqBkTHDhwE+EwxRYAhqm9z+UluIDYnWe9UrZ6cw2s\nr/Z1vc4rslF5swCAWnp2b04SHCGECmxEz2F6v2NdPT2Vu8z4Xx8Eq+jib7zWNb9AP84xdCAHwwTh\n1zNywkk0nETGSZR6cw3HW3JjICqsLWK82nazLQqRFX3NDEQMrrJkHKzmRb+tgmqXw28MTLnV9Vum\nPx+DNk/guP58r15UgTYhy6i3B0DHPAIrsXu0pcIUgOG4zQDgHT4obD8C0b6RW7pXUM/2pd7Ty/8b\ngmY3mwzo9jq9dtplacOxno2IGra6xWN5RTYshrILoNuejfSca5TH9Jz8UwDG0J3Dfwza0q1edbnL\ndNE9qf9v1K0AeI5SOksIqeghBUALhZjb7CJ2oW213Wc/scKSeRA0fKDO7wbKv3G2wSVM4CeE4Iab\n2ALadN2tzSGcybzVUe4luMv1QvBy3SYfbG0O4ULRcaa2L5EN05tV53cjU/ifnvXsSOVAGI1+h1ue\nHH6Putxl+kCYY56sOZTAQuyhg2bLexApKbdIdOZ3Ah3v9gV2EgU3IXEbEOrLNbC74byBoNttu2i8\nztWXa+DO5kbH8ihF1gsjpcupL8jZYBI7YhfapOGZluOROxlGgrtfsVUI9d6lQRdBkcLL06ZXqEC0\nyHohJyhIwiD20EE+N4/uDItkotzyk8Bt+A0hGMLiFLMFulf8Otuj3YqbwwleuAmjOdwQRJTtxNWw\n1cDtBySIyAYJGZh/ZP32hVSGDSSBid2jVQj/PPogsHR0R682q60wFdSr9RNCMODxbqtK5/KSfbnG\n8sMPQT1ft3Mbtrp5sUCwgS+h3qreF8yIChvINQ5WH7Ff0lpDnDcbdgdVr35k+f8ki60hVGMLdznW\nCyK4PLCK+9jCXZ6TELxENkhclsebBTr7wlpAm4JLmR6SbmIPHUjs8cpEcEr9MhjMVZEh3gtA2Akg\nS4iBpz0vWGwNK1zAUs4CizcrwwZrl9iFNpupAQi25TgvLGledjm1pO/Nznb0xWac8Jqa61XOIraA\nfdwWAG4UJh3jt25E4ekCnaGOG4VJ2zpBvFhAjIja3b1Y+4JE4kbsoYNC7mbcJjCjDP+I+5ig3tS6\n7KLvUMKb/a+tnMcj9hkVhh1WW8y2GkQhsrwhAwNzXxCZ0sUa/kpKyqOEjdiFdrEmdlaYyDit9QvU\nrnyuqw7LUnkibl1ZxNYqTO+58W+7z2USuqiEl+V8Zlvt3osVEesXBPF27fqCGzJssLaJPXTAyqBS\n8NykkQdRuy4A3iEEIHgYAWCbQeYVTug6r0X8eEIMrG2ywppREIUnC7D9iMoJChIWYhdaQlQA4ndg\nFUVHrDYTwipQJljFFnDeUddgIFeDmuEXvDjCC6y2snixkYksR19YFd6sCmTE+TlrjthDBz3F67Gd\nm/ULYHgt2dGvO7clIIRg1BERSgCAMxt/xHQbHgeGXYZtZzb+yLV+lCLLQnb068K9WZk/u3qJ/dJW\n60PC2wyrw7YufNy1XJTYstbzGii7c/J9y/9bhS0unGww22qGZTAQECuyLNdRufgnTG1JJEACQgeq\nGq8JrLHasrKEGw3vHwUR8Vreek7hhKKDHVahY43n+oFV2K228gx2RS2yANBmXAR+VYQNJIGJXWjT\nRAZsO0CKFlsAgQTXCycx5BFgUZ6yaIHlqccqsmVlCTNMNdnhuQuTqV3pI3ah7SneAMC2LQhP5gHP\nzrisXu3Azq+he4a7f3hElGdfMkOsXh0OtghOlGGGK5teTIXIGgzs/Jp3m6vImyUUyNTl9Fq/xB6j\nZVmPNinUZ8aYB0B4vrg8osEzmLO9uok5xhkHhm3rsosYXHDbr28Fns8gDJE1rn99Zoz5GIlECq0O\ni/dRn30XAPbcyTDElqfu4MLo8v9mUYsTJzvMttrB+yMTpsgCK33BsV0Ob1ZmG9ij7wt2kBCyR//r\nuAKVW90AZWP663sJIYf17W24bQMSEDoIE57wQViwxGsNeMIDPGEHK1aR443pBjkXL7zpWDz1/Yqs\nZ7urKGQQM0eNnWYJIa8AeBrAPh91/Za9BGCnvkfYegBHAdzPcFwXsQttIX8TwLq4zQDgHavt2fji\n8v9eGzl2tMsptgC7gLrVn1xvtzNzN25iyCLCIrxkq61+8l2jFFlzXwgCrze7VgbCjB1ojee62O3h\nreu3TGenaQvx5Xo8thnELrQAX4Cddyour1frJraEtDqe84otgFC8W6O+gXEcFbCoelShBsPWsAUW\n4B/4ssPaF5bblt4sAGzQvTyDcUrpuGNte0YBzFpemyaEjFm3CXer67eMUnrKJLIAcADAEz5sA5AA\noa03BuM2gZnFqw8jP/B6x2s8YguE691aj9sw9SAWei5zHRsHg5kljEzdj2v9v/V1LA+8IusUMrDr\nC6tZZIlKka0zeyw3KKUPBDzlekF1/ZYBWN52fC+AY6adb3lsA5AAoY0CkV6tHWGKLeBfcDNQbT3d\nJBB0Kqyf40WJrCjW4iAYIWQ/ALcUE0PQptG9maCTwLnV9VsGYHnb8ScJIfsJIccopQ9x2gYgAUKb\nzVYBsItOVNiJbX7wjGP9sMUW4Bfcat952+N52hCFlzBabfXbjh1+QgVeImvtC6vZmxUJRwihAhvx\ncrg1d6xLCIHPslEAeymlT+ovfxuAkXnAYxuABAhtPruAKIRWRAZCacNPXcv9iC3AHrc1YBXchfKv\nPduwI4gI+/VU3WwN0m4YIgt09gU/IuvHm10rA2FAhxACWL6FP255Pk0pnXWr67cMWhzWPOd+FMCs\n7uHCzTY7YhfapRrbnHEzotemdcLq1c6dO4B1d/yV6zG8Ygv4824Bb8G95a2P4tqO/893u1HiZGuU\nAguwhwtY+oIkMI8RQg4COAVgDMBjprJDAI4BGGeoy11GKT1OCCnroQ4AeAjAhxjb7CJ2oY0SP16t\nnwXC/YotwO/dAvGGBMIgqNCHLbId54rIm40dCijsg2FiTqndihu348ctZfs46vote870dNxS5nic\nHZ5Caxp1M5R73JL2YK47BuABaIHiBwE8Ybjazu23keSFv4EVsVVytm/bFj9iCwQTXKDTy23lgqcx\nRQXJz6VKYJXcbKQiu5bCBqsRFo+WaQaEPgXtASPYrSfwHoP7CKO+qMxmTrP9hw/8xmqHlBow6r2Q\niBnjC+xXcP2KLaAJbmvbt2EkzyXR0zULa3Xrd3y3EyQv1m9mwehtf+P7nJK1h+vvq90MCABOMyBG\nsZLQCwCvABj1mgO8VNvAZmkCqJ7/tK/j/H6Zh7ILgUSkdGllcWpjrQDzI0q8zm+2lZWgn4/v66LU\n8Na5A9zHpTJkIBGCl0fLPANCH8F7yPTSA9BG6VzvtymNPmzg16ttNdf53tQxqHdrwOPlKs1+13LR\nmQdBxNvLVjNBZ3YFyY81wgWtZjKmjUvSgZfQcs2AsMRjD8BhJE4fydsPAMPDQzhx+r8AALaO/AR9\npWs4O7EXAFDur2D3tu/iX898AQCQyTTw4F1fwenKJ7CwtBEAsHP0bzE7dy9uTL0bALBp4zHkcnO4\neElrY2DgLLZs+j7Ovv649oazC7hz91cwef6zWNJXyd8yehg3px/A/Ky2XsTQxu9BUZq4fuWjAIC+\nwd+gvOHHqNc2YuK1LyGbm0Fp59cxV/kM1KbmsA/uOozqjd9DY+4eAEDvphdAaRZLb38YAFAo/xKF\ndadAz38GbShAfgrZ7c+iVfk00NbWE8iM/i3UyQ+ALtwGAFA2fQ9o9UO9/n7tc1v3C6zvfwPtix9H\niypoF66jtuUF9Ex8EkTNAQAWd34dhWsfRHZpK5T6EJTqJmTqG5Cf1t5bY90v0S5dRunKwwCAdukq\napteRG/l0wAIAIrF0a9jZPL9yFQ3AQCqm19AproF+RltxarG+pNoF26gdFV7b62eS6iP/BC9lT/V\n+oHSxNKOZ1C8/DAy+mdc3fI8sgu7kJu7FwBQH3oZNLuA4jVtIJc0+wE1h96JT2ptZKpY2v4sSm99\nFEpDE7XCjm9Cnb0Prbl3aJ/P8A8BpQX17d/X2hh4DcrQCbTPa3YgN4fsjmfQmvgk0NSCKOt3/Xcs\nTn7Y8zrdPP8Z7ZoUrmNgxzcxe24/Mu0ezAPYdvtTaDbKmHjtSwCA4S3PodUcxPSk5mcMrH8ZfQOn\ncWXiswCAfPEK7hr9Bl597XGoqhZrveeuJ3Hx0h9jXr/WO7YdQbW2CdcmP6C1ectP0ddbQWXiUWRA\nMNB3CXfvOIoTZx4HKAEIxXvueQqvTuzDzYWtWpujRzC3sANvTb4XwMr3CZhAUIhKkakFn9K9ViGU\nOq81QAjZC+CAPhvCeG0GwIfcknN1IZ22jNrZ8jv3lej/eHEjn9UmgqZ58Xi2rWYfsqYBpqDblfvx\nbu1w8nJJqwSajX5XWz842SpiXQIg+Cwv68CXtS+4ETRkEHQgrGfzxMmgU2IH+rfQB+//PFPdH/74\nPwc+32rDqwtwz4DQB8EqLCILAI0AAz5RM3vj33U8DzobSNQUTyNWaRWl3Mz9DkckD7OtTu/HD2Vl\nSbjIAt19wYm4RVaSDFy7gVVQ7WZnWBfKhebJHtef7/UyoNUq8drcQdCOyPNFWNBvWc2IEFuRc+rN\nIpWbv01Yu2EylF1AaWFUmLgC4gTW6fra9QWJxAmW9C6m2Rm6CJ8EOqanVQAwebZxEnR6rvFlDBJK\nCDJY5kSWqF3CFSRlTBSixNQOYXcJAtYukFkGEgNPoWWdnaEPhBFwUsjPwceqYx2ImJLLIra3bH7e\ntdxvRoIZs1AEFV1l4w+6XnMTOZEizCumdrbyIPSugEFkvfqCCJGVYYPVQwKm4HJrc2h4ia2qj+y7\nIcK7NQjs5ap8lzdMT9MTTlsNohZYA7e+sCo9WQpk6jLrwC+xd4l6Y0BIO6J+/d2+JFNvf4S5HZHL\n5n/qvTEAAAsRSURBVBnxRl5RUSc/KMyGsOGx1e/n4Qbv9XLqC6JEVnqzq4sEeLTJQ9SmjiK9WwOR\noYU0EdYi3CJ/EFelJysRQuxCm8suAWCfFeSGyOUT7cS2v3zSV1thCC7gLbpk8DdCzxcmVlvD3t0g\nqMBa+4JIkZXe7OpjVQmtaKxiO7D+FefKDIQluIC9MM14LKadJMrrX0YmZHEFxHmw5r4gPVmJF7F3\nEdGLyoj2BgaVlS/S5Qr/QiJ2uOVnCuXCIx3xzLC9RFasNpWVpeUpr2Eh+jO/XDnQ0TdEkVRvlqgU\npNZieki6id2jDYMwdmAIw2sJ08N1gkVsg8R+kyLmQLj7eIWxFFJSRVYSnNiFVlFaSPrC3wY9heuh\ntGsWBJGim/Fpbxxi6ddWO8IUWOMHd1JwX5Ai2w3npgOOdQOUjUFL8i9D28rmkGnPsEMADkJb3fAV\naGvCOG5yELvQlgpT8LPwtxdheLW37/rq8v8ishLsECm6Azu+GdScyAhqaxShGPNdjbkvSEKDadMB\nhrp+y14CsJNSOksIWQ/gKABjUY5zlFLmSQAJiNHyb87Iimgv4bU3/tzUdviDIEZc0a+IzJ7b710p\nIfDaav5swhZZu2tt7gvB25ferBWeTQfc6vot09lp8qCnEYDYPVpKw1UrkZ5ty2aKqvEFDMvDNbCK\nCYu3SxOwrgErXrZGMnhowe2H1K4v+DvHqhTZDbp3aDBubHHFAfOmA251/ZZRSk9ZwhQH0LmDTFlf\nNGsWWljhL902OYhdaKMgiu3JoxJcAzvhiXJQLWziEFaDqNK1UiWylEKpNVhr3xCwHi3PAihudf2W\nAeiI4R4zViXUMcdyp6GFGRzXJY1daHuLkwA2hX4eEWJ71x1PMZxH+xuV4JqxitO63f8NilJPtAAb\nNhu2xgmPwLL0BfdzpUhkBaJvCuC2YashaNPQBqHMOAmjW12/ZQCWF8t6khCynxByzNgEwey96tt4\njRFCyk5ebexCW2+KWeuAhaBie/nqH2Dbre6rNq2cS/sbh+Aa3Hj7DzC85Xkm71C0GPN6pIatUePX\ne+XpC93nXJsiCwAcIQSeTQcc6+pLtvopGwWwl1L6pP7ytwEc1l8vA3jaGEQzHZfc0EGrHa23FURs\nb968CwDfl8v8RY5adJfm2e2N81Yd4LM1KCJCA376gnbutSuyPJiEEID9pgPQNhmYdavrtwxa/HbI\nZNIotM1mK/pmB4dNx+2Bx7rbsQttHBidPey4bfd5V/6P09Ndi8Q9TVYKrC+YNh1gqMtdRik9Tggp\n66EOQBvw+pBeNksIqZjKdsFhI1oD180Zo+Cd9w7QE8eGvCuGBI/Y3py/DQP9b4Zkh/g2l+ZvQ09I\n9oomDFvDFFeevhC3yIrYnHGwMEJ/d/Mnmep+f+IpuTmjhdg9WpXGOyuMx7tt6ltWh2NH53MRwtsK\n0V7RiLA1Sq+VpS/ELbCS5BD7hIVGMxkrdw0qBc8vxtW3H3ItF4mRJG9+8DI9GZ29QeG1VcTnEwS3\nvsDSlyRri9g92qQRV/yWBScxWc3x3rhjqzxIcZU4EbvQJnU9WjvB3TD0clzmuOIkRiNDLycizcwJ\ns91mW9OAuS9IgZV4EbvQZjNVJFFoDcxfotrg6Rgt4adssjfpIlZO2Wd767o30CsFVsJI7F+/aj2+\njANezlc+l6r425uVz8ZtAjNpsNW49oNKAb9589G4zYkWlQLVGttD0kXsHm1aMYttEuO5EjGk5UdV\nkmxiF1pFaSItC3/39bxt+7rdlzEJ4lsqXYnbBGaSYCuPqDr1BYnEjtiFtlSYRhgLf4fBvaPfYq6b\nBPHdtfMbkZ4vCFHbGtRT5ekLEknsMdql2nDcJjDz87PBFns2x/jMj7B49bXHQ2tbNGHZGtZnHrQv\nSNYWsXu0HLtBxE67nQ+lXZYvvh9vWFXTE1/0Y2uc8dOw+kJiUVVQOdDlm9iFVsLGah+UyYCs+vco\nWbvELrS9pWtIS4z239zz5bhN4OK99/w/UEg6xCttn23a7JXES+wx2lrDusB5cnnj4h/FbQIXabI3\nTbYC6bNXEi+xC227nQ6PCwBm50fjNoGLNNmbJluB9NkriZfYQwcSiST5UEqhysEw33gKrWkXSGMF\n8nGnvXF46hoUCzMANnCaHQ937XDdrSJxpMneNNkKpM/eNCJKe/yWWdo/TCk94Mc2gM2jPWpsQqbv\n1f40gH0C6gIAVDXHYEIyWKiOYLDvYtxmMJMme9NkK5A+e1OKKO3xWwb99T0A9gM4YHqZS+tcY7SE\nkDFoW/ICWN7lcU/QumYazT6vKonh0rX3x20CF2myN022AumzN22I0h6/Zaa2y3qdWdNr3FrnNRg2\naj6BzrR+oiB1JRKJxA1R2uO3zGCPzRbn3FrnFTro2vNcRF1990hjB8l6z2akZDHS/30DgBtxW8FO\nmuxNk61Ayuy9I2gD83T6xWPNI6yDKUX9dtpgnFI67ljbHlHa47fMCBkctynisQ2At9BOA7Amujqd\nhLmu/qGPA1p8Iy07ZqbJViBd9qbJViBd9lpEzxeU0j8QZMt+aNtzO3GMUnoc4rTHV5k+2DXtMMDF\nYxsAb6Gt2DVg40rz1pVIJGsQDs9WiPYQQuCzbC+A9YQQ48e0rP9IHOe0DYCH0JqMAbCs8sctz6cp\npbNedSUSiYQVUdoToKwjf09P7xo3PXe0zQ6W9K7HCCEHsZIv9pip7BCAY9DDAB51neCN3cRJmmwF\n0mVvmmwF0mVvmmw1I0p7/JYZWQf79f8PAniOUlrxOs4KoZSyvmmJRCKR+CD2tQ4kEolktSOFNsUQ\nQkYJIQcJIXv0v0xLoRFCDodtmyRaCCHHGOr46i+S4IQeOgh7rYQYbR0D8AC0NI8HATyhx24igxBy\n0jQNsAzgaUqp65RnPTfwGI14awvea6uP+i5jHZwIG599AdD6w3NR9QX9eo4COOx1Tf30F4kgKKWh\nPgCcNP1fhjZHOHDdOG3Vy/abnu8BcC5iW8egCab5tRmPY8r6ca71EtAPDgLYa6p7MkzbRNhreX44\nBnup6P4iH+IeoYYOolgrQRSc5x8F8ITp+SsARiO+FfMz5dluOmHocPaDMoC/oLoHS7X0nfsjMXTF\nBt6+eCAFt+FyinyMhB2jTdNaCczn18XqIdNLDwCYpRGGOcA5DdBlOmEU8FzbBwBUCCF7TbHEqFfZ\n5u2LhwCcJ4Ts15Pan3CoFyfc00Yl4gh74e9Q1koICa7z084Y3AGw5QyLhHkaoMd0wijg+WxHod3m\nHqeUzurTR0/CfdqmaHj7wrju0RrL6B1Ht1DHDfe0UYk4wvZoQ1krISR8nV/3YJ6lEQ/WgG8a4BiA\nB0weV1n/PypPkeezrQCoGD8K+t/RiL1arr5ACDlIKX1SD3EchpZInzTkFPkYCdujTdNaCdzn12/H\nK1RbBCNSKN8URdfphBHA2w+sRO0dMtur94FTpjrjhJBdhJCxuEVMTpFPDqF6tNaOZicGxiCCV92w\n4bFVfz4GrRMf1593pCNFxGNGXiS0VCTrFMWPmSsTQsr6tEFEGfvk7AcVALPGc/1vhUaYOsfZF6ah\n3TG4thEWhJAx0zU9pPcFA2sfcOsvkhCJIo92DNqIrd1+PUehpZyMe9WNAlZb9S/eOcvhFUpplHHE\nVMHZD0ahxTvPQYvNHo5SaH3YuxcrHnAZWnxZ3pJLlpFrHUgkEknIyCm4EolEEjJSaCUSiSRkpNBK\nJBJJyEihlUgkkpCRQiuRSCQhI4VWIpFIQkYKrUQikYSMFFqJRCIJGSm0EolEEjL/CwxhSM1Rlh+s\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4019511890>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUYAAAEACAYAAADP1t+BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztndtvXFeWn3+7LrxfihRJkZQluim3ux0nxlgtB2i4g3Zg\nCh60EQdtqD1IXvJkCUGCgTEDSOgAecjTQPoPRD/mIfBIsREBdiKImrgBGwbGkg0YceLx2GxdWiIp\nSlWkSIqXuuw8nH1KR8Wqc/au2zqLXB9QkKpq16mPp4qL+3bWUlprCIIgCE9IUAsIgiDEDQmMgiAI\nFUhgFARBqEACoyAIQgUSGAVBECqQwCgIglCBVWBUSl21aDOtlDqjlJox/2Ya1xMEQWg/Kmwfo1Jq\nBsA0gAtaaxV6IKVuaK1/Yf6fAfC+1vp3zZQVBEFoB6GBsdxIKR0WGJVSxwCc01qfCDyW01oPNUdT\nEAShfTRrjnEawErFY1kTMAVBEFjRrMA43KTjCIIgkJNq0nGyACoXW2oGS6XUKQCnAKCzq+MXhw6P\nAwCKyU1oVUSq0AcAKCXyKKbWkd55MiLPd+SQyg9A6aR3P/0IiVIHksWukGNsIL3j6WloFDpWKo6x\nikSxC8lSJwCgkHoMQCNV6DXH2EExuel5qBK0KqGQXkUqPwilE+VjJIvdSJQ6zDE2ACikCj2eV2Ib\npeQW0vlBz0MVUUg/QmonAwVlfrYVJAu9SJTS5hjrUDqJZLHb/GxbKCV2kM4PPHWMyvOTLPR5x9AJ\nFNKPdh1DJ/JI5fu9n00VUEyv7T5Gvh8J7X09Cuk1qFI64hw39jl5n8tqcz6n8jlu3eeU2hks32/4\nc6pxjGZ9Trf/X/aB1noUDXDiX3brh9mSVduvv9m5orX+80bej5pmBcZ5VAmEWuuvqjXWWs8CmAWA\nqX9yQP/lhz9tkkZr6dg8iJ3uJWoNazj5cnIFePm+98K1W40e42G2hM//16RV257JmyONvh81dQ+l\nzfacDLA7ACqlpgHMWR3H9AY4EPxrzQFOvpxcAX6+ghuhgVEpdUwpdcb8/5zZvuNzDsA7gfvv+vsY\nAZwE8K6NgD904MBA7iVqBSc4+XJyBfj5Cm6EDqVNT/ArAOerPPe7Gm0By96iIAhCHCG/JLCY3KJW\nsGYt8y21ghOcfDm5Avx8BTeatfhSNzqRp1awZrt7kVrBCU6+nFwBfr6NUoTGammbWqNtkPcY/e0I\nHBhZeJ1awQlOvpxcAX6+ghvkgVEQBCFukAfGkipQK1iz3XWfWsEJTr6cXAF+voIb5IGxmF6jVrAm\nO/4ptYITnHw5uQL8fAU3yAMjp42yE7d4ZVHj5MvJFeDnK7hBvirNivCUlPGDky8nV4Cfb4MUNZAt\n7Z+fWQKjCyo6d2Ws4OTLyRXg5xsjzCXDJ+FdEHIMwKzWujJtYWRbk9bwuGmaAXBJaz3fDEfywJjv\nyFErWLMwdZFawQlOvpxcAX6+MeNiINv/dQDvA6g1NxHWdkZrXb4qTyl1AcDpZgiSzzEmGe1jHF58\njVrBCU6+nFwBfr5xwfTysv590/ubqbPt6VbVliIPjH4uOQ50bo1RKzjByZeTK8DPN0a4ZPuPansO\nwB+VUqdMjtezzZLkE5UEQSCjAIUVkyDYghEz7PWZNTlYAbds/6FttdazpsfoD5/nsDuQ1gV5YCww\n2sf4YOIatYITnHw5uQL8fNvMA6318RrPuWT7D22rlDpj5hjPmx7jVQBH6/DdBflQWpnU7hzo3Byn\nVnCCky8nV4Cfb4xwyfZfs63J+/pV4LFZAJeaVYCPPDD69Sk40L/yIrWCE5x8ObkC/HzjQlS2f4fK\nAFl423dCj18v5ENpQRD2He+aygD+3sRgtv9z8IbEs2FtTa9x2gyhAW/I/UGzBMkDYzG5Sa1gzaOh\nb6gVnODky8kV4OcbJ8Ky/btUBtBaX2qVI3lg1KpIrWANp83oAC9fTq4AP99GKeokHhb7qDXaBvkc\no19XmAMHln5NreAEJ19OrgA/X8EN8sAoCIIQN8gDY4lRzZetnnvUCk5w8uXkCvDzFdwgD4zF1Dq1\ngjW50c+pFZzg5MvJFeDnK7hBHhglUW3r4OTLyRWw980Velts0h4KSCBb7LO67QXIA6Mg7FX8oLhX\nguN+QgKjAzrBp3AXwMuXkysQ7VsZDCU48oI8MHLaD7Z45ENqBSc4+XJyBcJ9awVBCY58IA+MqfwA\ntYI1IwtV82nGFk6+nFyB2r5RwU+CIw/IA6PSSWoFa9LbLqnk6OHky8kVqO07lNoIfV3U80I8IL8k\nUBD2GkOpjao9Q85BsagTWCn2UGu0DfLAmE8/olawZnnyCrWCE5x8ObkC0b6VwZFzUGw2TawSaH0c\nV8gDY6LUQa1gTffGFNY6+GRV4eTLyRWw8/WDowTFXTSrSqDLcZwgn2PklKi2b/Xn1ApOcPLl5ArY\n+0pQfJpmVQl0OU49kAdGQRD2Fc2qEuhyHGcih9KO8wHHAPhFcDIALmmt58OOzylR7erw19QKTnDy\n5eQK8PNtlIJOuGw1akeVwJZuY7CZY3QZx8+Yql0w7S/gSWnDqnBKVFtkNizi5MvJFeDn22baUSXQ\n5TjOhA6l6xjHn/YL2djCKVHt8P1fUSs4wcmXkyvAzzdGNKVKoONxnInqMdYcx9cQOAfgj0qps+b+\n2SptBEHYp5giVuX71aoEAshqrVfC2kYdp1GiAqNT11RrPWt6jP7weQ67AytMZa9TADAyNoSJm+8A\n8AoM5Tty5bTxWz33kBv9vJziSScKWDzyIUYWZspXHixPXkH3xlR5lXB1+GsUUxvlv+ibvXeweuA6\nxm//FgBQSm5h6fBljN57A6mdQQDA/UOfoPfR8+hdew4AsDLyJbQqYGj5lwCAx303sZb5BslCHyZu\nvoNiegP3D32MsbtvIpn35l2WnrmM/pWX0LP+LAAgN/oFlE4h8+AVAMBG/w/YGPgeY3d/AwAodKxi\nefIKDt55CwmzMr945CMMPjyO7o3DAIDs2GdIFnoxmH0ZALA++B02e29h9N4bAIB8ZxYPJuYwfvtt\nqJL3US5MXcTQ8qvoejyJZKEPHZsHkd4ZwkDuJQDAWuZbbHcvYmThdQDAdtd9ZMc/9c6xVoDSWJi6\niOHF19C5NQbAKy7fuTleLhnais9JqxJUKd2Uz+ngn94CgJZ+TlqVyt/bRj8nAHh48A8t+5xiSMNV\nAi2eawilta79pFInAZzWWp8IPJYD8Hq1HqNS6ow/x2iC31mt9dEwgSP/dFD/9aVX6vVvK6qUhmaU\ncZyTLydXgJfvey9cuxEy52fF5IsZfeoDuzo3/+WfXW74/aiJ2q5jPY5XSs3gSZlDmFWoS1HL5+kd\npylJUvzeDBc4+XJyBfj5NkpRJ7Ba6La67QVCA2NlAKw2HxBYbMnC686GHkMQBCHu2GzXsZoPMJOh\n02YIDXhL6R9EHVyj9lA+bpSSW9QKTnDy5eQK8PMV3IgMjKbH5/f65iqe+13F/UuuAoWOplzz3RaW\nDl+mVnCCky8nV4Cfr+AG+SWBnBLV+quMXODky8kV4OcruEEeGDklqvW3jXCBky8nV4Cfr+AGedox\nQRDiT1EnsJLfGyvONpD3GPPpVWoFa+4f+oRawQlOvpxcAX6+ghvkgTHBKB9j76PnqRWc4OTLyRXg\n5yu4QR4Yk6VOagVr/EvRuMDJl5MrwM9XcIM8MAr7CykfKnCAfPGlkHpMrWDNysiX1ApOxM3XD4rV\n6qDEzTUKbr6NUtIJrOf5THs1Sgx6jHyufNGqQK3gRJx8K3uKlffj5GoDN1/BDfLAmGI0tPLTW3Eh\nLr61hs9PlReNiast3Hw5Yi4xPqOUmjH/1sw4E9ZWKXVMKXXK3M6YnA+hkA+lhb1N1JyilBcVQmhW\nmVXnkivkPcZSYodawZrHfTepFZyIg29U0POfj4OrC9x8udGsMqsG55Ir5IGRU5XAtQyfgvBAfHxr\nBcfg47Vcc4Xe8i1OxOXcxpQRpdT1wO1U9Et20awyq8CTkiun/ATaUW9OHhjTeT7XnPop87kQJ9/K\n4Fh5P06uNnDzbZSiVniU77S6wVQJDNxmo45fhWaVWfWTZv8NvOHzaZtjyxyj0DaGUhsNzSnKfCR/\nTI8trNzJVa31HJpXZjVYcuW8ef+rEQ70gVGrErWCNcU0r1/KOPrWCmxxdA2Dm29ccOg9NqXMarWS\nK0qpoyGVTgHEYChdYJVE4mNqBSc4+dZyHUptPHWLC5zOLUdcyqpEtK2r5Ap5YEwxmmMcu/smtYIT\nnHw5uQL8fJnyrr83EcBJ7C6r8k5UWxMA54P7GGFRcoV8KK00eWy2xq9NzAVOvpxcAX6+jVLUCus7\n7U344lhWJaytc8kVPlFJEAShTZAHRk6Japee4VUAiZMvJ1eAn6/gBnlgTBb5pEvvX3mJWsEJTr6c\nXAF+voIb5IExUeqgVrCmZ/1ZagUnOPlycgX4+QpukAdGQRCEuEG+Kl2I0d60KHKjX1ArOMHJl5Mr\nwM+3UUpaYTOfptZoGzHoMSpqAWuUJv874gQnX06uAD9fwQ3ywJgq9FArWJN58Aq1ghOcfDm5Avx8\nBTfIA6MgCELcIA+MxcQ2tYI1G/0/UCs4wcmXkyvAz1dwgzwwlpJb1ArWbAx8T63gBCdfTq4AP1/B\nDfLAyClR7djd31ArOMHJl5MrwM+3UbRW2NxOW932ArK0JghCLDHpw07CSw5xDMCsqecS9pqrWusT\nVR4/GbwflViCPDBqVaRWsKbQwee6boCXLydXgJ8vU6yrBJp0Y9OoUjDLpBqb11pfMjkcrwGId2As\npB9RK1izPHmFWsGJRn2fqvvc4o34++3cCuFUq/xngl9VTDkEvzRq8DgZAL/XWg/5xwHwi6j3j5xj\ndCl6bdqfDN6ijp/acapqSMrBO7wKIDXT17VKn2tlv/18boWquFQJDOM4vES1JwMxbDrqRTY9Rpfu\nrHOXVTG68iVR7KJWcILKNxgQbQtYybmNN1or5PNJ2+YjJlb4zNZRKdClSmAY0/DmJ+dMr/M6gBto\npBiWS3e23i6rwAebIFetlyjV/fYdD7TWx6s90aIqgWHMw+usrQDlGDatlJrWWs/XelFUj7Fmd7ZK\nMZlyl9W85hiAS2FvDgD5jtBFplixeOQjagUnGvV1DWZhQ+eo4Ljfzu1+pUVVAqOOU0lk0IkKjC4R\n2rrLav5qnAKAkbEhTNz0ato8GvoG+Y4cDiz9GgCw1XMPudHPMXHLG7nrRAGLRz7EyMIM0tue2vLk\nFXRvTKFv9ecAgNXhr1FMbWD4/q8AAJu9d7B64DrGb/8WgLehfOnwZYzeewOpHW8P5f1Dn6D30fPo\nXXsOALAy8iW0KmBo+ZcAgMd9N7GW+QaTN/8tiql1FNMbuH/oY4zdfbNc+2PpmcvoX3mpnKcvN/oF\nlE6Vr6nd6P8BGwPfl/e/FTpWsTx5BQfvvFUeli0e+QiDD4+je+MwACA79hmShV4MZl8GAKwPfofN\n3lsYvfcGACDfmcWDiTmM334bquR9lAtTFzG0/Cq6Hk8iWejD/UMfI70zhIGcl1h1LfMttrsXMbLw\nOgBgu+s+suOfeudYK0BpLExdxPDia+jcGgMAPJi4hs7NcfSvvBj6Of3MFKEvqjx+OPQ/ceT+v0D3\nzhAA4ObYp+hfO1rzc9KqhKXD/6Mpn9NB36OFn9PBO/+6XK+o0c8JAB4e/EPLPieOmNKn5fvVqgQC\nyEZt39FazyulVpRSGROXMvB6kKEdNqW1rv2k1/s7HdwXpJTKAXi9SsnCGQAXtNZHA49pAEfDJI4+\nP6X/8vJPwxxjw8TNd7Dw7N9Sa1hD4Vur1xjV+5Rz2zree+HajVpDW1u6jh7Sz/zNv7dq++Nf/OeG\n3w8oT+X5daGf2seolLoIb9g9W9H2HIDzeDIk94PoaQA/wuuoXYgKjFE9Rtei15XwGScLTWEotbEr\nOMr8olAPdVYJPF/lOPMAzrq8d2hgdOnO1ttlLaTWXXxJyY59Rq3gBJVvMDjaBkU5t/FGlxQK2+Tb\nntuGzU/6rtmG43dnK4teXwXgT6j+DsDvlVJ+l7Xqtp4gSltvASAn6biXjxpKX9deopxbIU5EBkbH\n7qxzl5VTlcDB7Mt4PPCP1BrWcPLl5Arw8xXc2D99Y0FgSrUFLZm3bS3kaceKjPIxrg9+R63ghI2v\n66V+rWIvnluBL+Q9xlJih1rBms3eW9QKTkT5xiUoAnvv3O45NKDz5P2otkH+k6bzA9QK1vgbdrkQ\n5hunoAjsrXMr8Ie8xyi0n7gFRSEcmU9sP+SBkVOi2nxnNrpRjKjlW7kJOw4bsvfKuRX2BuRDaU6J\nah9MzEU3ihFhvrWCX71BsdFe6F46t0Gkd84T8sCYNkkGODB++21qBSeifCuDYKM9xUaCwF47t8CT\n8yHBkR/kQ2lO+JlRuGDj6w+rmzF8buQYe+3cVgZD9jkptYLaJu9HtY3985MKNWlGT5H1L32TqdVD\nlJ6jGy5lVZRSx5RSp0y7i7XKF1TWhKkF+Z9pTvniFqYuUis4wcmXkytQ2zcq+MkfESesyqqYgHk8\nkIJsBl4Oh6MV7Wbg5YE9HfXG5IExWeijVrBmaPlV5ByyqlCv9tr61vplvbk1gr9fnirf7+vYBgCs\n73Q+1e7EhP1VILXeq5prnINIrXNbLe1a5fNCNI5VAqfh5Wjwk9lcBzDtZ/oyx8uY41mlQiQPjIlS\nmlrBGj/bMhdsfIMLBMFf2ptbI1jJd2Mp14/CQy9zte4sAcBTc02Hnluu+nqX96rmGtY2DoSd21rB\nMY4/R4yxLqti0iOeCDx0HMBKRXbvGVOkz+rNyQOjQEe1BQIAWC10YyXfjfV8FwrbKXQtJpFeA0od\nT1LEJXaAjWdL2Mynn3q9/8vfyJByLyxcVAZHbv670EAib13Rs+1VAivyvp5GID2i6Wk67QcjD4yc\nEtU+PPgHagUnwnxrBa5gULy3MVC+Prb7oUbno1K53cpPqufRdFlgCAY83zVs4SJOwcXmu9DMFX9m\nkFUJNMf/QGt9ydy3qg1TCXlg5JSoNr0zhJ3uJev21L8QYb7VhnvBoPgo34nNfBpqO4GkN7WIrqUt\npBayWP8zbxhZSuvy6wZTm3U5+oEjvTOEpXT4H8k4BRnb70JcfONCK6sEmp7hvF/rxXAMwLBSyg/U\nGRM858KqC5Bv1+GUqNav4MaFKN+h1Eb5F3cotYHB1CYy6U30pbcwkN7GaO+TQJXa8oJgabgfyS2v\n5+gwtAp18F2jgkicggy37wI3qhTb21VWJbh9x1+sCRTAOmmOc0lrPevfzGOzjRbDEmKEax0VV/zj\nDqY2sZr2/mD9Q984AK9Xv3G4G4D3eMc6sLmtyu1tjgsgct5NFi6EAFZlVUzQvAEAgcWVeQCX/Dsm\niJ4y/z+DiJr35IGRU6Latcy3ZO9dmfTBJlC4+PrH/Emnt8qMTmCl2IPuvm3k+zqxdih6cFEtANpe\nkx105bBwQfld2C/YllUxAS50+GLmGM+jShXBapAHRp3IUytYs929SPK+1XpQNsHRxjcYwDLJxxhO\nruNA0htCPyz24chwDrdfBB7lupF+6H1dSp3am1/sK5SPU+26a5eeX6Vr3BcuqL4LVKgSkNhufOqE\nC+RzjKl8P7WCNSMLr7f9PcNWeaNWgKN8awXFTGIb06ktHE0/xH88/Hf4i59+hX/+8z/i4J8tAVOP\nkXkuh0PPLeOVo7fw68kfmpKpp5prXIMiQPNdENoHeY9RoCEsqA4nNAYTnQC2sZJcx1BqA8/05LyV\n6r40RnvXMdmzikx6E892PWiftCC0CfLAWFKF6EYxYbvrftvfM+wSs8ihdIhv8LjVjj+d2kK25A2d\nMsnHAICB9DbQi5YERYpz2wjcfAU3yANjMb1GrWBNdvxTq3bNXj2uFhxtjm3rC9TuQT4seteyD6Y2\ngR7gT4+HWtJTdHGNA9x8BTfI5xg5JaqduLUrsccuwkoGNEIwENoG3CjfyuPkCr1YKfYgW+zDQ3PL\nFvuwUuwp73N8pifXkuGzzbmNE9x8BTfIe4ys0OGrcvWuHtvifJwIX/+YYcF8pdjzVNuWLYhYuMYK\nbr6NovfXqrQERheUrvlU1OoxyQpriG+Qam7BgNgWLF1jAzdfwQnywCiJaltHo742UwH1BPyqx91n\n51aIN+RzjElG+xiHF1+r+VxYgKDajxfm2wya+XO12rXZcPMV3CAPjAlN3mm1pnNrLPT5Wtf+UhHl\nGyc4uQL8fAU3+EQlJgQXM+J85YYgOKEBRlfvNgx5YCww2sf4YOKaVbu4BERb33qpd1Gp2mta7dps\nuPkKbkQOpV1KGFa8zqpMoWJU86Vzc5xawYlGfdsZ4PfbuRWiqaN86oxS6qRS6kKwfKptadUgNnOM\nF7XW500CyFl4JQyjfiC/TGEkyWKXTbNY0L/yIrWCE83wbdei0n48t0IkLrHnGoDrpqTBDQAXgadL\nq2qtzwO4AC+PYyihgbFaCUMAtUoY+q9xKlMoxJ+4LSoJe586Ys9PAnVdsoHH/dKqPuXSqmHvH9Vj\nrFnCMOQ1M2F1GSopJuurFULBo6FvqBWcaKZvPZckurCfz+0eZEQpdT1wsxo9VuAUeyqKXZ2GCYYm\nFkWVVt1F1OKLUwnDesoUalV0aU4Kp83oQPN9W9lL3O/nNu4ojXJRNAtqVgl0wCn2AOW6MCfxpNIg\ngPDSqrWICozWJQxdyhSavyCnAGBk9AAmbr4DwPsrnO/I4cDSrwEAWz33kBv9vHzBvk4UsHjkQ4ws\nzCC97WksT15B98YU+lZ/DgBYHf4axdQGhu//CgCw2XsHqweuY/z2bwEApeQWlg5fxui9N5DaGQQA\n3D/0CXofPY/etecAACsjX0KrAoaWfwkAeNx3E2uZbzB+57fId2RRTG/g/qGPMXb3TSTz3tacpWcu\no3/lJfSsPwsAyI1+AaVTyDx4BQCw0f8DNga+x9jd3wAACh2rWJ68goN33kLCzLMuHvkIgw+Po3vj\nsHfyxz5DstCLwezLAID1we+w2XsLo/feAADkO7N4MDGH8dtvQ5W8j3Jh6iKGll9F1+NJpHeGsXj4\nI6R3hsrFm9Yy32K7e7GcaHW76z6y459651grQGksTF3E8OJr5b16DyauoXNzvDyv1orPKVnow93p\n/9qUz+ngn94CgJZ+TmN330TRlP5t9HMCvHKsrfqc4kQry6eaAHjeLLRc1VoHe4q7SquGempd+5pP\n0219X2v9i8BjOa31rpQ4pipXUPwCvOgcWqbw6PNT+i8v/zTKMxZM3HwHC8/+LbWGNWM3/w3+4ZnL\nVm2p5wy5nVtOvu+9cO1Goz247onDevrf/ZVV2/977q8afj/H2DMN4KRZXPHXOXIAjvqxx4xmUVFa\ntSahPUat9VeBqltVSxjC9BIro7BS6oJNDdkSo12jWz33qBWcWO+yr0tCXV+F27nl5ssNl9gDbz7y\nQODl0/DmEf2g6JdW/crcPxnVa7TZ4G1VwjAg7FSm0B+OcCA3+jm1ghN3D3zp1J4yOHI7t9x8mWIV\ne7TWc0qpTGCR5wSA14FyAA0trVqNyMBoW8Iw8JhTmUJuiWq5DJ8A4Gd3/5X1UNqHKjhyO7fcfBum\nBCR22vuWLrGnogc4G3g8srRqNcgvCRRaRxLFqkEuKp0Y9bBaEKiRwOiATvAp3AXU9g0rsOU/D9iX\nZmhGEG3VuXWpbe0Ct++C4AZ52rG4bScIY/HIh9QKToT5NqMWdDPZS+dW4A95YEzlB6gVrBlZCL0a\nMnZE+VYGQcrh8147twJvyAOj0klqBWv8zcpcsPH1gyH1nOJePLcCX2SOUSAPikL8URpItnlVmhLy\nwJhPP6JWsGZ58gq1ghON+rYzYLbq3LbqZ+D2XRDcIB9KJ0od1ArWdG9MUSs4wcmXkyvAz1dwgzww\nckpU6ydA4AInX06uAD9fwQ3ywCgIghA3yOcYOSWqXR3+mlrBCU6+nFwBfr6Nst8WX8h7jJwS1RaZ\nrd5y8uXkCvDzFdwgD4ypQh+1gjV+UlUucPLl5Arw8+VIsyqU1nMc8qG0IAhCDS76iWqVUtfhVQnc\nldErSKBC6elGjkPeY+SUqHaz9w61ghOcfDm5Avx8udGsCqX1HAeIQWDkNFezeuA6tYITnHw5uQL8\nfNtM26sEGqpVKK3nOPSBMb1jNW0QC/xCTVzg5MvJFeDn2zAaSG5rqxtMlcDALbLESRWaVaG0rova\nZY5REIS20YoqgREVSp2rDQIxCIwatasUxo1ScotawQlOvpxcAX6+ccGh9ziPKgGsylAZ8OrBDCul\n/MqEfv2XOcfjlCEPjIWOyDLUsWHpsFv9FGo4+XJyBfj5cqOZFUrDjlML8jlGTolq/QLqXODky8kV\n4OfLlHf9/YcATmJ3lcB3go1NpcAz5v9nTBCMOk5VyHuMnBLVpnYGqRWc4OTLyRXg58uRZlUoDTtO\nLcgDoyAI8UeVgNQWn/WARiEfSufTq9QK1tw/9Am1ghOcfDm5Avx8BTfIA2OCUT7G3kfPUys4EUff\nXKG3aknTOLqGwc1XcIM8MCZLndQK1vSuPUet4ETcfIMBsTI4xs01Cm6+ghvkgVHYH1TrJVZ7TBDi\nAPniSyH1mFrBmpWRL6kVnIiLb1gAzBV6MZTaiI2rLdx8G8VLVCuLL22Ez8nWqkCt4AQnX06uAD9f\nwQ3ywJhiNJwaWv4ltYITcfENK2HqPxcXV1u4+QpukAdGYX9QLTi2s261ILhAHhhLCT4Vdh733aRW\ncCJuvsFAWBkU4+YaBTdfwQ3yxRdOVQLXMt9QKzgRR99avcQ4uobBzVdwg7zHmM7zueb04J/eolZw\ngpMvJ1eAn2/DlDSSWyWr216AvMcoCIJQDZMd5yS8BBDHAMzWSEZb+boLWuvTgfvHAJRzNQK4pLWe\nDztGZGB0kQsIZAC8AuBslIBWfP7CFNO8Fgs4+XJyBfj5MqVZVQJntNbnA20uVDy/C5seo5WcqdBV\nru9gBK8iPI05CqySSHxMreAEJ19OrgA/X25Uq+5nYkrYa3ZVCTScVkpZ9TZ9QucYHUsPTgM4G7h/\nHcB0VHHs2oZrAAALQUlEQVTrFKM5xrG7b1IrOMHJl5MrwM+XIc2qEgh4SW3/qJQ6ZUoenK3S5imi\neow15SoFTCryE4GHjgNYiYrSSpOv/1iTzPPZjA7w8uXkCvDzbRSlgeR20bb5iBld+szWUSmwWVUC\nobWeNR00f/g8h91x7SmiAqOTXMV84mnUSCFuovYpABgZPYCJm16G8kdD3yDfkcOBpV8DALZ67iE3\n+jkmbnkjd50oYPHIhxhZmEF621NbnryC7o0p9K3+HACwOvw1iqkNDN//FQCvMPrqgevlcpel5BaW\nDl/G6L03ylmY7x/6BL2Pni9nTFkZ+RJaFcpXNzzuu4m1zDdI7wxj4uY7KKY3cP/Qxxi7+2b5F2Tp\nmcvoX3kJPevPAgByo19A6RQyD14BAGz0/4CNge8xdvc3AIBCxyqWJ6/g4J23yqnXFo98hMGHx9G9\ncRgAkB37DMlCLwazLwMA1ge/w2bvrXJa/XxnFg8m5jB++22okvdRLkxdxNDyq+h6PIn0zjA6Ng8i\nvTOEgdxLAIC1zLfY7l7EyMLrAIDtrvvIjn/qnWOtAKWxMHURw4uvoXNrDADwYOIaOjfH0b/yYss+\np2ShD6qUbsrn5K8Yt/JzShb6yt/bRj8nAHh48A8t+5wIeKC1Pl7tCYIqgVBKnTFzjOfN+0dO8Smt\na1+rrJQ6CeC01vpE4LEcgNfDqmyZN89WFqmpxpEXM/qv/3vVcxg7EoUulFJ8qsNx8uXkCvDyfe+F\nazdqBSpbBgae0a8c/w9Wbf/uf/+nht/PDJnf99c3zGM5rfVQlbYn8XTQ9BdX5uCNemGCrd/+HIAP\nwmJY1DjWufSg6dLO2wRFAEgWu22axYL+lZeoFZzg5MvJFeDny43KGFOtSqC/fqG1vqS1nvVv5rFZ\nM4LNwttNE3r8SkIDo4ucuX8MXk9xztw/GXZ8AEiUOqKaxAZ/+MUFTr6cXAF+vkxpuEqgiWHz/sKL\nef6DqDe22a7zrjmYv4+xUu4qgFkTNG8YKf/5eQBWPUdBEIQgTawS6ByDIgOjrZzptio4UmCUYSU3\n+gW1ghOcfDm5Avx8G6akkdiyXpVmTwz2yjjHUjKU5nUFJSdfTq4AP1/BDfLAmCr0UCtY42/p4AIn\nX06uAD9fwQ3ywCgIghA3yANjMbFNrWDNRv8P1ApOcPLl5Arw8xXcIJ8oKSV5bJIFgI2B76kVnODk\ny8kV4OfbKEprJLbz1Bptg7zHyClRrX+ZGBdcfHOF3qdurlS+3vUYcm6FOEEeGAVBEOIGeWDUis/e\nqEIHn9yRAC9fTq4AP1/BDfLAWEg/olawZnnyCrWCE5x8ObkC/HwFN8gDY2onNI9trDh4h1cBJE6+\nnFwBfr6CG+Sr0orRlS9+Pj4uuPjWKmvartfLuY05JY3EpqxKC4Ig7FvIe4z5Duv6NOQsHvmIWsEJ\nTr6cXAF+vhxxrFB6DsAZeCULrsNLsD0feP6pFIhRGXfIe4xJRvu6Bh/yyDTuw8mXkyvAz5cpF7XW\n501+11l4FUpr8aPWWmmth7TWJyqC4hmgHAznAPw+6o3JA2OilKZWsMavxcIFTr6cXAF+vtxwrFAa\ndpwMgN/7PUSt9UqwXEItyAOjIAh7jhGl1PXA7VQdx3Atn5pRSp1USs0opc4FKgsch5fB23/ujBmi\nh0I+x1hIrVMrWJMd+4xawQlOvpxcAX6+DVMqAZvWeQ1qVgl0wKlCKQLzj0qpLIBrAH4BL8AeAzCn\ntV4xZV1vIKJKIHlgVDpJrWANp/lQgJcvJ1eAn29caEX5VKA81Pb//5VS6pjpNc7DK8634rcztaqm\nK8o9PwV5YORUJXAw+zIeD/wjtYY1nHw5uQL8fOOCX8XPAusKpdVKrZq2K0qpasEvciuMzDEKghA7\nHCuUzsOrJe0/NwNThM/0Clf8tn4vMqy3CMSgx1hklI9xffA7agUnOPlycgX4+TLFqkKp3zMMLPIc\nrWj7OwC/V0r9aJ7bVWGwEvLAWErsUCtYs9l7i1rBCU6+nFwBfr4NU9LAVnuz7buUT/Vr2dc4zjyA\nsy7vTT6UTucHqBWsGb33BrWCE5x8ObkC/HwFN8gDoyAIQtwgD4ycEtXmO7PRjWIEJ19OrgA/X8EN\n8sDIKVHtg4ma0xixhJMvJ1eAn6/gBnlgTO8MUStYM377bWoFJzj5cnIF+PkKbpCvSnNClXidLk6+\nnFwBfr4NUypBP96ktmgb5D1GQRCEuEEeGPMdOWoFaxamLlIrOMHJl5MrwM9XcIM8MCYLfdQK1gwt\nv0qt4AQnX06uAD9fwQ3ywMgpUW3X40lqBSc4+XJyBfj5Cm7ssxlkQRDqQWuNkn0+RvZEBkbHgjTW\nbX04Jap9ePAP1ApOcPLl5Arw8+WIazyxKXillLqgtT4d9d42PcaLfp4zk/32fdTOTuHS1hNllKg2\nvTOEne4lag1rOPlycgX4+TLFOp6YLDzzWutLJrXYNZjUY4E2MwBOAYgMjKFzjC4FaeotXsMpUe1A\n7iVqBSc4+XJyBfj5csMx9kQWvDJtsrBIUgtEL764FKRxLV4jCIJQC5d4YlPwaqZa9u9aRA2lXQrS\nWLc1CSX9pJLb771w+/84vA8h10YAPKC2sIeTLydXgJnvzxo9wJrOXrma/28jls27zNDXZ9ahpIGP\nS+wJLXhlhtBOF7dHBUaXgjTWbc1JmgW8uYMmVBRrC5xcAV6+nFwBXr4VQaoutNZ/3iSXVhTDqlnw\nyjyfjVoEriQqMFoXpHFsKwjCPqQVxbBM20r8QHgMwLBSyv8jljHBea7uKoGmDGH5frWCNDDROKqt\nIAiCLY6xZ14ptaKUypjeYrDg1VPBz2zXiQzONtt1rArSWLSthevcAyWcXAFevpxcAV6+nFyDuMSe\n0IJXJlieMv8/A+BSWI9Raa2b9UMIgiDsCcivlRYEQYgbEhgZY1bezgT2blWu4tV63YXoVgInlFJX\nLdrU9X3Zj7R8KN3qa60JXY/B21iaAfAKgLNhcxatQCl1I3DJVAbA+5X1dqu8ZgbelggV1q7ZtOK6\n11ZS53cB8L4PofNXTfacgbeP70LUZ1rP92XforVu6Q3AjcD/M/Cuf2y4LaWree5U4P4MgB/b7HoM\nXoALPpaLeE3GvC60XQy+B2cAnAy0vdFKt2b4Vty/QOCrm/192c+3lg6l23GtdbNwfP9pAGcD968D\nmG7z0KSeSzCdLotqFs2+7rXV1PFdPM1gWCqX7DrQ6jlGTtdaW7+/CS4nAg8dB7Ci2zjsh9slU3Vd\nFtVEmn3da6tx/S6eA/BHpdQps3n4bI12lDh9X/Y7rU5U25JrrVuE0/vrp+eQTsNuz2Yzsb5kKrgZ\ntuVW1Wnada9twvW7MGt6jH46qzlYZnFpIy6X2O17Wt1jbMm11i2irvc3PYQPdJsXB+B2ydQxAMcD\nPZqM+X+7emINXfcKb5qinb1Gp++CUuqM1vq8GfJfgLfxOG7IJbsOtLrHyOlaa+f3N8PTee1d9N5W\ntNslU5UJO60ui2oizbrutV1Y+5rvwFeBNrNKqaNKqWPUQUcu2a2flvYYK78Y1X55/UnrqLatxsXV\n3D8G70s3Z+4/tb2kTbzr70uDt7Wk8pKpd4KNlVIZczkU2jl35/g9mAfgX+/qL8bM6zZuhXL8LmTh\n9chDj9EqlFLHAp/pOfNd8Kn8DoR9X4QA7djHeAzeit6u/WBKqYvwthDMRrVtB7au5hflx4qXz2ut\n2zkPxgrH78E0vPk6/7rXC+0MjHX4nsSTHmYG3vyoDFEZI9dKC4IgVCCXBAqCIFQggVEQBKECCYyC\nIAgVSGAUBEGoQAKjIAhCBRIYBUEQKpDAKAiCUIEERkEQhAokMAqCIFTw/wFuDu+E63QqZAAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f401cad1110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl = plot(u)\n",
    "plt.colorbar(pl)\n",
    "plt.show()\n",
    "\n",
    "pl = plot(local_project(p, V))\n",
    "plt.colorbar(pl)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "RuntimeError",
     "evalue": "\n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to solve nonlinear system with NewtonSolver.\n*** Reason:  Newton solver did not converge because maximum number of iterations reached.\n*** Where:   This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2017.2.0\n*** Git changeset:  unknown\n*** -------------------------------------------------------------------------\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-14-c6af918a3cf7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m     \u001b[0mobstacle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0mh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterpolate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobstacle\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0msolve_Uzawa\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mpl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-12-d2d4361d090f>\u001b[0m in \u001b[0;36msolve_Uzawa\u001b[0;34m(tol)\u001b[0m\n\u001b[1;32m     43\u001b[0m         \u001b[0mtold\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     44\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 45\u001b[0;31m         \u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mform\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     46\u001b[0m         \u001b[0;31m#gap.assign(project(h-u, P))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m         \u001b[0;31m#p.vector()[:] = np.minimum(0*p_old.vector().get_local(), rho.vector().get_local() +\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python2.7/dist-packages/dolfin/fem/solving.pyc\u001b[0m in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    298\u001b[0m     \u001b[0;31m# tolerance)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    299\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mufl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclasses\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mEquation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 300\u001b[0;31m         \u001b[0m_solve_varproblem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    301\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    302\u001b[0m     \u001b[0;31m# Default case, just call the wrapped C++ solve function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python2.7/dist-packages/dolfin/fem/solving.pyc\u001b[0m in \u001b[0;36m_solve_varproblem\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0msolver\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mNonlinearVariationalSolver\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    348\u001b[0m         \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparameters\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msolver_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 349\u001b[0;31m         \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    350\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    351\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: \n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to solve nonlinear system with NewtonSolver.\n*** Reason:  Newton solver did not converge because maximum number of iterations reached.\n*** Where:   This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2017.2.0\n*** Git changeset:  unknown\n*** -------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "for d in np.linspace(1e-4, 1e-3, 11):\n",
    "    obstacle.y = d\n",
    "    h.interpolate(obstacle)\n",
    "    solve_Uzawa()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<dolfin.cpp.la.Vector; proxy of <Swig Object of type 'std::shared_ptr< dolfin::Vector > *' at 0x7f401e7570c0> >"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "assemble(form)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
